Exact solution of a stochastic protein dynamics model with delayed degradation 
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We study a stochastic model of protein dynamics that explicitly includes delay in the degradation. 
We rigorously derive the master equation for the processes and solve it exactly. We show that the 
equations for the mean values obtained differ from others intuitively proposed and that oscillatory 
behavior is not possible in this system. We discuss the calculation of correlation functions in 
stochastic systems with delay, stressing the differences with Markovian processes. The exact results 
, allow to clarify the interplay between stochasticity and delay. 
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I. INTRODUCTION 



Due to the small number of some molecules involved and to the uncontrolled environment, biochemical processes 
inside a cell usually need to be described by stochastic models P, 0|. Some important basic processes (such as 
transcription, translation or specific degradation) are indeed compound multistage reactions involving a large number 
of steps of similar duration, and, in principle, due to the central limit theorem, the time to complete such processes 
should follow nearly a Gaussian distribution, rather than exponential, with a well defined characteristic delay time [H . 
A description including delay is then needed to obtain a reduced model for this kind of processes. It is well known that 
delay can change qualitatively the dynamical behavior, allowing, for example, the appearance of oscillations in the 
evolution of the number of molecules and there has been a great interest in delay- induced oscillations in biological 
systems. 

Stochastic processes that include delay are analytically difficult due to their non-Markovian character. Most the- 
i , oretical studies consider a Langevin approach (stochastic differential equations)^, H| or systems in discrete timeQ 
(where delay can be accounted for by increasing the number of variables). None of these approaches is completely 
suitable to describe chemical reactions inside a cell since the former neglects the inherently discrete nature of the 
molecule levels and the latter considers an arbitrary discretization of time. In some cases, the discreteness can be a 
major source of fluctuations Q. 

In this work we develop a rigorous derivation of the stochastic description of a protein dynamics model that includes 
delay in the degradation, and solve it exactly. We find that the exact solution for the probabilities leads to equations 
for the mean values that do not comply with simple intuitive arguments and that oscillatory behavior does not exist 
(while it is usually believed to be present in this type of system). This clarifies and warns about the derivation of 
dynamical equations describing the evolution of the concentrations in cases in which delay plays a role. The exact 
\^ • solution is specially valuable for small system sizes, where approximated schemes typically fail. Due to the low number 
of molecules inside cells, this regime may be biologically relevant. Our solution extends the results of reference 
in which the authors consider a particular case which allows for a Markovian description in terms of suitably defined 
,_H ' dynamical variables, while such a Markovian reduction has not been achieved in the more complete case studied here. 
t-H . The paper is organized as follows: in the next section we define the model and in section Mil we derive the 
corresponding master equation. In section IIVI we solve the master equation and derive the probabilities and the 
(macroscopic) equations for the mean values in the case of constant rates. Some technical details are left for the 
Appendix IVIII In section [V] we discuss the details of the derivation of the master equation for the conditional 
probabilities and derive, again in the case of constant rates, the correlation functions of the process. Finally, in 
section fVTl we end with a brief discussion of the results. 



II. MODEL 



Proteins usually degrade through complex proteolytic pathways that involve several different steps (such as tagging, 
binding of auxiliary proteins, recognition by protease and destruction) [1J| . As discussed before, it is then natural to 
consider a time delay in this process. As a simple model for the dynamics of the number of molecules of a protein X 
we consider: 

C 7 P 

0^A, A^0, A^0, (1) 
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where, for simplicity and to isolate the effect of delay in degradation, transcription and translation steps have been 
lumped into a single stochastic process that occurs at a rate C. The delayed degradation (indicated by a double 
arrow) is modeled as a reaction that is initiated at a rate (3 and completed at a time r after initiated (giving rise to 
the destruction of the protein) . Besides delayed degradation, we consider also instantaneous degradation at a rate 7 
(which can take into account processes such as non-selective degradation, proteins going to the membrane or outside 
the cell, dilution due to cell growth, etc.). This model of protein dynamics was proposed in [l2j and it was thought 
that it can lead to periodic oscillations, although the analysis of [£j, in a particular limit snowed otherwise. To 
completely describe the process, one has to specify if a protein that initiates delayed-degradation at time t (and thus 
will disappear at t + r) can also disappear before the completion of this reaction, through instantaneous degradation 
at a rate 7', not necessarily equal to 7. The process is equivalent to the following two- variable system: 

c 7 P 7' 

->■ X A , X A 0, X A X It Xj =► 0, Xj 0, (2) 

T 

where we have split the proteins into two types: Xj are "infected" molecules that will die precisely at a time r after 
being infected (if still present) and X A are non- infected particles (so X = X A U Xi). In principle, we allow the rates 
to depend on n A , the number of X A particles, but not on rij, the number of Xj particles which are considered to be 



"inert" . One has to include an extra variable because delayed degradation is a "consuming reaction" [1J| , since as soon 
as one particle initiates this reaction the number of particles that can initiate it decreases by one, and so the state 
of the system changes both when the reaction is initiated and when it is finished, after a time r. In a recent work 
[Tol | , the stochastic description was studied in the particular case that infected proteins can not undergo instantaneous 
degradation, 7' = 0. In turns out that this particular case allows for a Markovian description in terms of suitably 
defined dynamical variables. Our treatment allows us to pose and solve a non-Markovian problem and it can be of 
interest to other processes not allowable to a Markovian description. 



III. MASTER EQUATION 

To derive the master equation of the process, we start with the following identity, valid for any stochastic process 
(with a numerable set of states, otherwise the sum should be replaced by an integral): 

P(n,t + At) = ^2P(n,t + At;n',t), (3) 

n' 

where P(n, t) is the probability of being at state n at time t and P(n, t; n', t') the joint probability of being at state n at 
time t and at state n' at time t'. In most cases of interest, P(n, t+ At; n', t) is O(At ) or O(At) only for a small number 
of n! (usually depending on n), whereas it is o(At) for all the rest. Writing explicitly the terms O(At ) and O(At), 
dividing by At and taking the limit At — > one can obtain the master equation of the process (differential equation for 
P(n, t)). For Markovian processes one can write the joint probabilities as P(n, t + At; n', t) — P(n', t)P(n,t + At\n' , t) 
and the conditional probabilities can be readily derived from the rules of the process, obtaining in this way a closed 
equation for the one-time probability P(n,t). For non-Markovian ones, in general, the equation for the one-time 
probability will depend on higher order joint probability densities P(n, t; n%, t%; ...;n rn ,t m ). In some cases (as in the 
one considered here) it is possible to explicitly write the joint probability densities that appear in the master equation 
and obtain a closed equation for P(n,t). 

In our process, n has two components (n A and ni) and the only terms O(At ) and O(At) that appear in ([3]) are 
associated to the following elementary processes schematized in ©: 

(1) Birth of an X A particle from the reservoir: (n A — !,«-/) — > (n A ,ni), with probability CAt. The contribution to 
the right-hand-side (rhs) of Eq.© is P(n A - l,n I ,t)CAt. 

(2) Death of an X A particle: (n A + 1, nj) — > (n A , Uj), with probability j(n A + l)At. The contribution to the rhs of 
Eq.© is P(n A + 1, nj, t)j(n A + l)At. 

(3) Infection of an X A particle: (n A + 1, Hj — 1) — > (n A ,ni), with probability f3{n A + l)At. The contribution to the 
rhs of Eq.© is P(n A + 1, m - 1, t)/3(n A + I) At. 

(4) Death of an Xj particle: (n A , m + I) — > {n A , m). This might occur by two different reasons: 

(4a) Death by instantaneous decay at rate 7' with probability j'(nj + l)At. The contribution to the rhs of Eq.© 
is P{n Al ni + M)7'(nj + l)At. 

(4b) Death after time r after infection. This is the only non-Markovian contribution. In order to account for this 
case, and according to the previous discussion, we need to write the corresponding contribution to the rhs of Eq.© 
as: 

^2 P(n A) ni,t + At;n A ,ni + 1,*;<S;Z„^), (4) 
n ' A ' n 'i 
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where 1 n ' A ,n' — Wa — l,n'j + l,t — t + At; n A ,n'j,t — r} denotes the event in which there were n' A and n'j particles 
at time t — t and one particle got infected during the time interval (t — t; t — t + At), and S denotes that the particle 
infected during that time interval survived up to time t. We use now the series of conditional probabilities: 

P(riA,ni,t + At;riA,ni + l,t;S;l n > Atn > i ) = P(riA,ni,t + At\nA,ni + l,t;5;X n ^ jn ^) (5) 

xP(riA,rii + l,£|<S;2^ in j) 

X.P{S\I n ' A y i ) 

xP(l* A ,n0- (6) 

Now, it is P(riA,ni,t + At\n A ,ni + 1, t; S;X v j lTl ^) = 1, as the particle that was infected in (t — r,t — r + At) and 
survived up to i, dies with probability 1 during (t,t + At), according to the rules of the process. The conditional 
probability of the second line of the rhs is the probability that there are n A , ni + 1 particles at time t given that there 
were n' A , n\ particles at time t — r and one Xa particle was infected in (t — r, t — r + At) and survived up to time 
t. Since the presence of the infected particle does not influence the dynamics of births, deaths and infection of other 
particles (remember that we have assumed that Xj are "inert" particles such that the different rates do not depend 
on the number of Xi particles), this conditional probability can be simply obtained by considering the process in 
which the dynamics of the infected particle is decoupled from the rest of the particles, i.e. as if removing that particle 
from the dynamical process. This leads to: 

P{riA,ni + l,t|5;2^^ )n j) = P{n A ,ni,t\n A - l,Uj,t - t + At). (7) 

Furthermore, as all Xj present at t — t have died before t (except the one infected during (t — r,t — r + At)), 
this conditional probability is, in fact, independent of n'j and, for translational invariance, depends only on the 
time difference r, so we can write it as P(nA,ni,t\n' A — l,t — r). The next term, P(S\L n i n i ), corresponds to 

the survival probability during a time interval of length t, or e -7 T . Finally, P(I„^ n ^) is the infection probability 
(3n' A AtP(n' A , n\,t ~ r). Collecting all the terms, the contribution to the rhs of Eq.Q is J2 n > A P(nA, nj,t\n' A — l,t — 
t + At)0n' A P(n' A , t - T) e -f' T At. 

(5) None of the above processes occur in the interval (t, t + At): (n A , rij) — > (ua, ni). 

P (nA,ni,t + At;riA,ni,t) = P(riA,ni,t + At\riA,ni,t)P(nA,ni,t) (8) 



1 - CAt - ( 7 + /3)n A At - jmAt - P{n A ,n I ,t+ At-S;!^^ \n A ,nr,t) 



P(nA,ni,t) 



= [I - CAt - (j + P)n A At - j'mAt] P(n A ,n!,t) - ^ P(n A , m, t + At; n A , ni, t; S;!^^). 

n A' n 'l 

Treating the last term in the same way as in the previous case (4b), one obtains that the contribution to the rhs of 
Eq.© is: 

P(n A ,ni,t + At;n A ,n I ,t) = [l-CAt-{~f + [i)n A At-in I Ai\P{n A ,n I ,t) 

- P{nA, nj - 1, t\n' A -l,t- T)(3n' A P(n' A , t - r)e~ 7 ' T At. (9) 

n A 

After adding up all contributions to the rhs of Eq.©, dividing by At and taking the limit At — > we obtain the 
master equation of the process: 

dP{nA,ni,t) = {E -i _ 1)Cp{n ^ nut) + {EA _ l)inA p {nA ^ lit) 

+{E A EJ 1 - l)Pn A P(n A , ni ,t) + (Ej - l)j'niP(n A , nj, t) 

OO 

+(Ej P ( n A> n i -l,t\n' A -l,t- T)(3n' A P(n' A , t - r)e^' T , (10) 

where E A j are step operators, E A jf{n A j) = f( n A.i + k) and P(ua, t) — y^ n P(n A , nj, t) is the marginal probability. 
As written, this is an equation for the one-time probability, and has to be supplemented with the appropriate initial 
conditions. The evolution equation for the conditional probability P(jia, ni, t\n A , n'j, to) will be considered in section 
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IVl For Markovian processes this equation would be the same as Eq. tfTOf conditioning all appearing probabilities to 
(n^,rij,<o), but for non-Markovian processes this is not necessarily the case. 

The next step is to calculate the conditional probability appearing in Eq. (fT0|) . Let P*(n A , m, t\n A , nj, 0) be the 
probability that there are n A particles of X A type and nj particles of Xi type at time t given the initial condition 
that there were n A and n® particles at time t = and under the condition that particles can not die through 
delayed-degradation. The key point is to note that: 

P{n A ,ni,t\n' A ,t-T) = P*{n A , ni ,T\n' A ,0,0). (11) 

This is so because, as reminded before, all Af-particles present at any time die before a interval of length r passes 
and their presence does not influence the dynamics of births, deaths and infections of other particles. Note that 
this equality is only valid as written, as in general P(n Al ni,t\n' A ,t — t') ^ P*(n A ,ni,t'\n' A ,0, 0), the equality only 
occurring for t' = r. For t' < r not all the infected particles present at t — t' have died at t; for t' > r some of 
the particles infected after t — t' may have died trough delay-degradation before t. Therefore, we need to compute 
P*(ua, nj, t\n' A — 1,0,0) (we simplify notation to P*(n A ,ni,t\n' A — 1)). It follows the same master equation (p~0|) 
without the last term. For the sake of completeness, we write down the final system of equations: 



dPjjiAVniji) 
dt 



{E A X - l)CP(n A ,ni,t) + (E A - l)Tn A P(n A ,ni,t) 
+{E A EJ 1 - l)/3n A P(n A , ni ,t) + {Ej - l)y ni P{n A , ni ,t) 

OO 

+(Ej - 1) ]T P*(n A , nj - 1, r\n' A - l)f3n' A P(n' A , t - r)e^' T , (12) 
n' A =0 

dP*(n A ,ni t\n' A -l) = , _ 1)cp .( ^ _ 1} + {Ea _ lhnA p* {nA , ni ,t\n' A - 1) 

at 

HEaEJ 1 - l)Pn A P\n A ,ni,t\ri A - 1) 

+{E I -l)~f'niP*(n A ,n I ,t\n' A -l) (13) 

The process defined by Eq. (|13[) is a standard Markovian process for which exact and approximate schemes 
have been developed to find its solution. Once we have obtained P*(n A , nj, t\ti' a — 1) by solving Eq. (fT5|) with 
the appropriate initial condition P*(n A , nj, 0\n' A — 1) = S nA>n ' i8 ni fl, we can replace in Eq. (fT2"]) and proceed to 
solve for P(nA,nj ,t). A very convenient procedure is to derive a differential equation for the generating function 
G(sA,si,t) = J2n A m S A A S T P( n A, ni ,t) in terms of the corresponding generating function of the process without 
delay-degradation, G*(s A , Si,t), see the Appendix for details. Note that the active variables X A follow an independent 
dynamics as can be seen by the rules of the process (the presence of Xj molecules does not alter at all the dynamics 
of X A ) or, from a more formal point of view, by noticing that, as derived from Eq.(|12p. the marginal probabilities for 
the numbers n A , P(n A ,t) follow a closed master equation of a birth-death process: 

dP(n A ,t) _ j 



dt 



= (E A L - l)CP(n A , t) + (E A ~ l) inA P(n A , t). (14) 



IV. SOLUTION IN THE CASE OF CONSTANT RATES 

Up to now, we have been rather general. The only assumption needed is that the rates can depend only on the 
number of active particles n A . We now provide the exact solution in the case of constant rates. In this case P* 
corresponds to a simple birth-death Poisson process. As shown in the Appendix, it turns out that if the initial 
condition P(n A ,ni,0) follows a Poisson distribution (this includes the case which starts with no molecules at all at 
t = 0) then the generating function at arbitrary time is G(s A , si,t) = e :E A(*)( s A-i)+2;7(*)(s/-i)_ This shows that the 
joint probability of n A and nj is the product of independent Poisson distributions, with mean values and variances 
(n A (t)) = <Jn A {t) — x A (t), (nj(t)) = of lf (t) = xi(t). It follows that the total number of X particles, n — n A + nj, also 
obeys at all times a Poisson distribution with parameter x(t) = x A (t) + xi(t). If we start with an initial condition 
different from Poisson-distributed, the time-dependent solution is not Poissonian, but this form is recovered in the 
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steady state with mean values (n A j) st . The mean values satisfy differential equations including delay terms: 

dx A {i) 



dt 
dxi{t) 
dt 



C — ax A (t), 
-j' Xl (t) + P(x A {t) - e^' T x A (t - r)), 



(15) 



with a = ft + 7. The solution with initial condition x A (t < 0) = 0, Xi(t = 0) = is 

-(l-e-*), 



X A {t) 

xi{t) 



(16) 



a— 7' 
CB 



l-e' 



(l-e T < 



< t < T, 
t > T. 



The steady-state values are (n A ) st = C/a, (ni) s t = C/3(l — e~ 7 T )/(a-f'). 

It is important to note that it is not possible to write a closed equation for the evolution of the total number of 
particles x(t) = x A (t)+Xj(t). In the case 7 = 7' we arrive at x{t) — C — jx(t) — (3e~~t T x A (t — T), and in the case 7' = 
we get x(t) = C — ^x A (t) — f3x A (t — t) which differs in both cases from the closed result x(t) = C — jx (t) — f3x(t — t) 
used in [12[ ■ A different process that leads to this closed form of the rate equations corresponds to delayed production 
and linear negative feedback: 



C-/9n 7 

=>■ X, X — > 0. 
r 



(17) 



However, in this case the creation rate C — f3n may become negative, and so the process is ill-defined If 



V. TIME CORRELATIONS 



We now turn to the calculation of the time correlations. For this, we will write a master equation for the conditional 
probability P{n Al m, t\n A , rij, to) and derive evolution equations for the conditional averages (n A ,i,t\n A ,nj,to)- As 
commented before, for a Markovian process, P(n,t) and P(n, t\n°, to) satisfy identical master equations but with 
different initial conditions. For non-Markovian processes, however, in principle one has to specify the whole history 
before to. In our system, the problem is that the dynamics from to depends on when the rij particles were created. 
We are mainly interested in the steady state. The term that causes difficulties is the one corresponding to delayed- 
degradation. To obtain the proper contribution to the master equation, and in the line of the previous arguments, we 
take a starting point similar to Eq.©: 

P(n A , nj, t + At; n A ,rii + 1, t; S;l\n° A , n T ,to) = 1 x P(n A ,ri[ + 1, t\S;l; n° A , nj,t ) 

xP(S;l\n A ,n°j,to), (18) 

with I — (J n , n , I n i >n i the event in which a particle was infected during the interval (i — t, t — r + At), independently 
of the number of particles present at that time. 

One has to distinguish the intervals t — to < r an d t — to > t. In the former, following the reasoning of Eq.([7]) we 
see that the first conditional probability in the rhs of Eq. p^)) equals P(n A , m, t\n A , n° — l,to), since, as commented 
before, one of the X\ particles present at time to will stay until time t, and its presence does not influence dynamics 
of creations and deaths of other particles. Next, we use Bayes' theorem: 

P(S;l\n° A ,nlt )= £ P&I*,^ K,n?,t ) = £ P«, T&tolfrZn^) pifo^o'"/ • ( 19 ) 

n A , ni n A , ni 

An argument similar to the one used in deriving Eq.([7]) shows that Pin^n^to^^X^ <n i) = P(n^,?7.j — l,£o|i^ — 
1, n'j, t — t) which can be written as P(n' A — 1, n'j, t — t\vP Ai n® — 1, to)P(n A , — 1, to)/P(n' A — 1, n'j, t — r). We use 
also that P(S;I n > ) = n' A P(n' A ,n'j,t — r)e -7 T /3At to obtain: 



P(S;l\n A ,n j,to) - V P(n' A - l,n^ - r|n° A ,n? - l,i )n^e^ T Ai 



f^'v-A ■i-a.-j '.^..a.^ ^P(n A -l,<,t-r) 

,i A ,n.j 

P(n A ,n°j-l,t ) 
P«,n?,£ ) 



(20) 
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We know that in the steady state, if the rates are constant, n A and nj follow independent Poisson distributions. This 
allows us to compute the ratios of probabilities in this expression: 



P(n' A ,n'j,t - t) (n A )s 



P(n' A — 1, n'j, t — r) n A 
P«,n°-Mo) _ n°j 



P(n A ,n%t ) (nr)* 

Which leads to: 



(21) 
(22) 



P(S;2K,n?,*b) = ^PK-Ui,t-r|<,n?-M„)^ V Ai^» (23) 

, , \ n l)st 

n A' n I 

= ^ e -7'r Ai M £in ! (24) 

where we have used n , P(n' A — l,n'j,t — T\n A ,n l j — l,to) = !• Putting all the pieces together we find that the 
delay degradation term in the master equation for the conditional probability for t — to < ^ 1S: 

(Ej - l)P(n A ,m - MK.n? - l,t ) el ? T _ l n° I . (25) 
On the other hand, for t — to > r, expression (TTg|) equals: 

^ ™/, ~ 1) t - r ; n A) n /> *o)e~ 7 T n' A f3AtP(n' A , n'j, t - r\n A , n j, t ), (26) 

which is the same as the delay-degradation term in the master equation Eq. (fT0)) but conditioning all probabilities to 
n A , mPj at time to , as happens in Markovian processes. 

Using (1231) and (|21)|) in the corresponding time intervals, we obtain the evolution of the conditional averages in the 
steady state: 

d{n A ,t|n^,n?,to> „ , . Q . . . 

= C - a(n A ,t\n A ,nJ,t ). (27) 

at 

d{nj,t|n^,n°,t ) _ (—y{m, t\n% n° v to) + P(n A , t\n% n?, to) - ^—^ n v t <t<t +r, 
dt \-j{ ni ,t\n A ,n%t ) + /3({n A ,t\n%n° I ,t ) - e-^' T (n A ,t - rln^n^to)), t>t + r. 

The steady state correlations, Kuv(t) = ( n u{t' + t)nv(t')) s t — (nu)st( n v)st can be obtained integrating the previ- 
ous equations and averaging over initial conditions in the steady state. The result is: 

K AA (t) = {n A ) st e- a \ (28) 
Kn{t) = !^ e ~i-:-^' T > °^^ T ' (29) 




KiA{t) = { " R Ja wl_^ ' _ t .7 - (30) 

a — 7 7 

K AI (t) = 0. (31) 

The autocorrelation function for the total number of particles n = n A +nj is K(t) = K AA (t)+Kjj(t)+Kf A (t)+K A j(t). 
Note that in the case 7' = the correlation for the X] particles decays linearly in time instead of the typical exponential 
decay (as can be seen taking the limit 7' — > in Eq. ([29)) ). For arbitrary 7' the correlation function for the Xi particles 
drops strictly to zero at t = r. 

An analysis of the previous expressions shows that, in all cases, the autocorrelation functions K AA (t), Kn{t), and 
K{t) decrease monotonically in time, so there is no signature of stochastic oscillations, as found in the particular 
case, 7' = 0, considered in [10j. These exact theoretical expressions are in excellent agreement with the results of 
extensive numerical simulations of the stochastic process (J5|) using a conveniently modified form of the Gillespie's 
algorithm [lg, 1 171 ]. This agreement provides an independent check of the correctness of the theoretical calculations. 
An important ingredient of the process is that particles which are infected at, say, time t and are then scheduled to 
die at time t + r, can nevertheless die instantaneously at a rate 7' during the interval (t,t + r). If this is the case, 
in the stochastic simulation one has to remove that particle from the list of scheduled events to happen at t + r. 
Otherwise, one is removing a particle that was already removed. Stochastic oscillations do appear if one, erroneously, 
removes twice these particles. 
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VI. DISCUSSION 



In this work, we have studied a stochastic model of protein level dynamics that includes delay in the degradation. 
The exact solution shows that no oscillations are present in the system, contrary to previous results, and in agreement 
with a recent analysis in a simplified version of the stochastic model that allows for a Markovian reduction [loj . 
This implies that the presence of delay in degradation alone cannot give rise to oscillations. We have also analyzed 
the derivation of master equations and the time correlations in systems with delay, pointing out the differences with 
Markovian processes. Our approach allows one to deal with more general systems that the one used in [101 ] and may 
be of interest to study other non-Markovian processes. Exact results are specially valuable in small systems where 
approximated schemes typically fail. They are also important to study the validity of assumptions and approximations 
to be used in more complicated systems. 

In a broader context, due to the important role of stochasticity and delay in many biological, physical or technological 
systems, the present work seems relevant as a case where exact results can be obtained, allowing to clarify the combined 
effect of stochasticity and delay. 



VII. APPENDIX: CALCULATION OF THE PROBABILITIES 



Let us define the generating function of P*(n A , tij, t\n A , rij, 0) as G*(sA,si,t) 



oo oo 



s A A s T j'P*(n A , nj 7 t\n A , n° l7 0). Using standard techniques, it follows from the master equation (fT"3 

ua=0 rij—Q 

that it obeys the partial differential equation: 



BG* BG* BC* 

— = [ 7 (1 - s A ) + P(sj - s A )]— + 7 '(1 -bj)-^ + C(s A 1)CT 



(32) 



with initial condition G* (s A , sj , 0) = s^s^ 1 . By the method of the characteristics, one can find that the solution in 
the case of the initial condition n A = n' A — 1, = (as needed in the master equation Eq . (pT0|) ) can be written as: 



G*(s A ,s I ,t) = [l + $(s A ,s I ,t)] n ' A - 1 eiq> \c [ dt' $(s A , Sl ,t') 



(33) 



where 



(s A , Sl , t) = (s A - l)e~ at + (sj - l)-^-f (e-^' 4 - e"" 4 ) . 

a — 7 V ' 



(34) 



It follows from Eq. ffT2l) that the generating function for the original delayed process G(s a > s i > t) 

oo oo 

J2 £ s'/sfP^m^) obeys: 



riA— rij—0 



^ = [ 1 (1-s a )+0( Si -s a )]^-+ 1 '(1-s i )^- + C(sa-1)G 
+ /3e~ 7T (l -sj) ^ G*(s A ,si,T)n' A P{n' A ,t - t), 

n'=0 



(35) 



We assume n A (t) = nj(t) = for t < so that both the master equation (p~0|) and the equation for the generating 
function (|35l) are valid for t > 0. The number of X A particles follows a Poisson distribution at all times P(n A ,t) = 
EA&LJL e - x A(t)(s A -i) with x A (t) the average value given in Eq. p^|) (assuming the initial condition is n A (0) = 0). 
Using Eq. (|33|) . it is possible to perform the sum over the n' A variable to obtain the equation: 



BG . , . , BG .. BG 

- = [7 (1 - SA ) + ^- s ^_ + y ( i- Sl) _ 

+ C( SA -l)G + (3e-^' T (l~ Sl )x A (t-T) 



(36) 



x exp 



x A (t - t)<$>(s a , Si ,t) +C dt' <f>(s A ,si,t') 
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One can check by direct substitution that the Poisson distribution G{sa, sj, t) = e x A{t){sA-V)+xi{t){si-i) j g a solution 
of this equation if x_a(i) and xj(t) obey the differential equations (|15[) . Therefore, and n/ follow independent 
Poisson distributions at all times (assuming 71^(0) = n/(0) =0). It follows that the total number of X particles, 
ri = ha + ni, also obeys at all times a Poisson distribution with parameter x(t) — XA{t) + xi{t). If we start with 
an initial condition different from n.^(0) = 71/(0) = (or ua Poisson-distributed), the time-dependent solution is not 
Poissonian, but this form is recovered in the steady state. 
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